The goal of study 1 was to…
library(tidyverse)
library(ggdist)
library(ggside)
library(easystats)
library(patchwork)
library(brms)
df <- read.csv("data/study1.csv") |>
mutate(
Date = as.Date(Date),
Participant = fct_reorder(Participant, Date),
Screen_Refresh = as.character(Screen_Refresh),
Illusion_Side = as.factor(Illusion_Side),
Block = as.factor(Block),
Education = fct_relevel(Education, "Master", "Bachelor", "High School", "Other")
)
# plot(estimate_density(dfraw[dfraw$Participant == "60684f29dbfe1bb2059e5e27_rkqoy", "RT"]))
# Dear participant, thank you for participating in our study. Unfortunately, our system detected multiple issues in your data (such as implausibly short responses, and random-like pattern of answers), which make the data unusable for us. We understand that you might have been in a hurry, or had some other issues, and kindly ask you to return your participation. We hope to open-up more slots in the future would you be interested to participate again.
outliers <- c(
# Half of the trials have very short RT
# Prolific Status: REJECTED
"60684f29dbfe1bb2059e5e27_rkqoy",
# Block n2 with very short RTs
# Prolific Status: RETURNED
"5e860198a846e30497df5189_6e43s",
# Error rate of 46.2% and short RTs
# Prolific Status: RETURNED
"615f319bb341cf2f306d858d_qsaq5",
# Error rate of 46.2%
# Prolific Status: RETURNED
"614f36fd81c78b7a125c4262_6ax4g",
# Error rate of 47.9%
# Prolific Status: RETURN REQUESTED
"61280140171ec546e87ed8cb_qdlgy",
# Error rate of 42.1% and very large RT SD
# Prolific Status: RETURN REQUESTED
"5d398380b37ab1000111fac3_2nxxh"
)
We removed 6 participants upon inspection of the average error rage (when close to 50%, suggesting random answers) and/or when the reaction time distribution was implausibly fast.
For each block, we computed the error rate and, if more than 50%, we discarded the whole block (as it likely indicates that instructions got mixed up, for instance participants were selecting the smaller instead of the bigger circle).
dfsub <- df |>
group_by(Participant) |>
summarize(
# n = n(),
Error = sum(Error) / n(),
RT_Mean = mean(RT),
RT_SD = sd(RT),
) |>
ungroup() |>
arrange(desc(Error))
knitr::kable(dfsub) |>
kableExtra::row_spec(which(dfsub$Participant %in% outliers), background = "#EF9A9A")
| Participant | Error | RT_Mean | RT_SD |
|---|---|---|---|
| 61280140171ec546e87ed8cb_qdlgy | 0.479 | 262 | 296 |
| 615f319bb341cf2f306d858d_qsaq5 | 0.465 | 263 | 372 |
| 614f36fd81c78b7a125c4262_6ax4g | 0.462 | 630 | 611 |
| 5d398380b37ab1000111fac3_2nxxh | 0.421 | 507 | 1679 |
| 5e860198a846e30497df5189_6e43s | 0.402 | 492 | 725 |
| 6104696d120a50b0ec51aa1f_m56p5 | 0.300 | 723 | 579 |
| 61572ca3e91309ebe876a9db_8cqnp | 0.287 | 659 | 333 |
| 5d9091ff391a60058a7711b5_dvz9e | 0.269 | 578 | 172 |
| 5bb511c6689fc5000149c703_r4kjk | 0.262 | 716 | 271 |
| 5ab308940527ba0001c15a9f_rnclh | 0.262 | 588 | 206 |
| 6106b7157977b80c497314f8_d7ukm | 0.260 | 718 | 1294 |
| 61088fcf393bde58359957b3_ol8jp | 0.255 | 837 | 711 |
| 60684f29dbfe1bb2059e5e27_rkqoy | 0.251 | 599 | 1326 |
| 611eb7284490ba01cddfbe98_om6zf | 0.246 | 699 | 414 |
| 5bc84a5f0760100001b3b9de_4cxmv | 0.246 | 560 | 197 |
| 60d129f2a122e84175a56425_z2w8h | 0.243 | 693 | 232 |
| 5d7389f193a945001a3ea504_nhua6 | 0.238 | 1160 | 1660 |
| 60dae077e62179eb469e32a4_b4pte | 0.227 | 748 | 243 |
| 5c6b0a27ffc824000191c7d8_5ajt1 | 0.225 | 780 | 427 |
| 5f2f0ac775b91d2e1865c1cc_xkcs9 | 0.219 | 785 | 836 |
| 5ff46a1a99e7cfb2994f7958_f2zg0 | 0.216 | 506 | 150 |
| 5f19559b9665f700090276c4_hsmss | 0.215 | 738 | 375 |
| 5c8ab0f10de08f00016e43a1_pyvgt | 0.213 | 1076 | 557 |
| 5ecd37ee75736a068808fa6c_7fgo5 | 0.210 | 870 | 489 |
| 6166a03f5063db088c458b73_m7w8f | 0.207 | 804 | 378 |
| 606cd013f538ed55e02069b5_vr3v7 | 0.206 | 652 | 367 |
| 5f480e566265722a9b2b2660_0bola | 0.205 | 511 | 147 |
| 605b60879326739b05897042_bpsyp | 0.203 | 627 | 223 |
| 609193e5a0cea97bf00ac6e2_a6zcr | 0.202 | 1133 | 982 |
| 610b0a1bf2434edb31592209_3f1wq | 0.202 | 869 | 424 |
| 5f08583a3d61a604d606c517_o75t7 | 0.201 | 720 | 298 |
| 55eab7fd748092000daa98f2_f10fa | 0.198 | 1110 | 738 |
| 612ba4de7e2b127155f4bb03_ph1f8 | 0.196 | 867 | 652 |
| 5e04595a4fa02aefdb9c9ced_n3rey | 0.189 | 983 | 830 |
| 61114f10ae21c59c0ed3d106_jw6v8 | 0.187 | 711 | 195 |
| 5f14886922a7d20725a22cde_9awyt | 0.186 | 803 | 397 |
| 60a6dd8779e3de1097e5d50a_t4wyc | 0.185 | 846 | 765 |
| 5c73e5d89b46930001ee7edc_ydo84 | 0.182 | 1045 | 1082 |
| 5e84f2663a34f20c3907e237_rt0oo | 0.181 | 1001 | 562 |
| 5ebde9baaefecd1325ef23c7_lphsv | 0.176 | 1307 | 1091 |
| 5d59a9d909f4300001de0c3b_l125y | 0.175 | 1146 | 900 |
| 563bb259be9cac0005aab7ab_te1z4 | 0.174 | 703 | 243 |
| 5fd97d5b40332e276ea58209_xmckd | 0.173 | 1046 | 918 |
| 60ba6031b6dde7c5b869bf74_gqplc | 0.173 | 616 | 381 |
| 5dfae1f373d7248254527108_0rb1e | 0.171 | 927 | 550 |
| 60b8e0ec34553723e3d6504d_ju18r | 0.170 | 769 | 312 |
| 5d5051e31025380015dc59b8_dwrdh | 0.157 | 848 | 364 |
| 60366cfe9748fc2b0ccbc9d0_ox8hj | 0.156 | 712 | 383 |
| 5bce155e561901000121006f_49472 | 0.142 | 1109 | 863 |
| 5ccc3dd7a758ba00133c0763_lwl1g | 0.139 | 895 | 816 |
| 5a0b46e0844c7a00015e4d13_jedw6 | 0.124 | 741 | 331 |
| 5eb0205cac7ad4085dc32a50_5xekt | 0.092 | 884 | 509 |
temp <- df |>
group_by(Participant, Illusion_Type, Block) |>
summarize(ErrorRate_per_block = sum(Error) / n()) |>
ungroup() |>
arrange(desc(ErrorRate_per_block))
temp2 <- temp |>
filter(ErrorRate_per_block >= 0.5) |>
group_by(Illusion_Type, Block) |>
summarize(n = n()) |>
arrange(desc(n), Illusion_Type) |>
ungroup() |>
mutate(n_trials = cumsum(n * 56),
p_trials = n_trials / nrow(df))
# knitr::kable(temp2)
p1 <- temp |>
estimate_density(at = c("Illusion_Type", "Block")) |>
ggplot(aes(x = x, y = y)) +
geom_line(aes(color = Illusion_Type, linetype = Block)) +
geom_vline(xintercept = 0.5, linetype = "dashed") +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
labs(y = "Distribution", x = "Error Rate") +
theme_modern()
p2 <- temp2 |>
mutate(Block = fct_rev(Block)) |>
ggplot(aes(x = Illusion_Type, y = p_trials)) +
geom_bar(stat="identity", aes(fill = Block)) +
scale_x_discrete(expand = c(0, 0)) +
scale_y_continuous(labels = scales::percent, expand = c(0, 0)) +
labs(y = "Percentage of Trials Removed", x = "Illusion Type") +
theme_modern() +
theme(axis.text.x = element_text(angle=45, hjust = 1))
p1 | p2
# Drop
df <- df |>
group_by(Participant, Illusion_Type, Block) |>
mutate(ErrorRate_per_block = sum(Error) / n()) |>
ungroup() |>
filter(ErrorRate_per_block < 0.5) |>
select(-ErrorRate_per_block)
rm(temp, temp2)
# RT distribution
p <- estimate_density(df, select = "RT", at = c("Participant", "Block")) |>
group_by(Participant) |>
normalize(select = "y") |>
ungroup() |>
mutate(color = ifelse(Participant %in% outliers, "red", "blue")) |>
ggplot(aes(x = x, y = y)) +
geom_area(data = normalize(estimate_density(df, select = "RT"), select = "y"), alpha = 0.2) +
geom_line(aes(color = color, group = interaction(Participant, Block), linetype = Block)) +
geom_vline(xintercept = 2500, linetype = "dashed", color = "red") +
scale_color_manual(values=c("red"="red", "blue"="blue"), guide = "none") +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
coord_cartesian(xlim = c(0, 3000)) +
theme_modern() +
theme(axis.text.y = element_blank()) +
facet_wrap(~Participant) +
labs(y = "", x = "Reaction Time (ms)")
p
# ggsave("figures/outliers.png", p, width=20, height=15)
# Filter out
df <- filter(df, !Participant %in% outliers)
p1 <- estimate_density(df, select = "RT", at = "Participant") |>
group_by(Participant) |>
normalize(select = "y") |>
ungroup() |>
ggplot(aes(x = x, y = y)) +
geom_area(data = normalize(estimate_density(df, select = "RT"), select = "y"), alpha = 0.2) +
geom_line(aes(color = Participant, group = Participant)) +
geom_vline(xintercept = c(150, 3000), linetype = "dashed", color = "red") +
scale_color_material_d("rainbow", guide = "none") +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
coord_cartesian(xlim = c(0, 3500)) +
theme_modern() +
theme(axis.text.y = element_blank()) +
# facet_wrap(~Participant) +
labs(y = "", x = "Reaction Time (ms)")
df$Outlier <- df$RT < 150 | df$RT > 3000
p2 <- df |>
group_by(Participant) |>
summarize(Outlier = sum(Outlier) / n()) |>
mutate(Participant = fct_reorder(Participant, Outlier)) |>
ggplot(aes(x = Participant, y = Outlier)) +
geom_bar(stat = "identity", aes(fill = Participant)) +
scale_fill_material_d("rainbow", guide = "none") +
scale_x_discrete(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0), labels = scales::percent) +
see::theme_modern() +
theme(axis.text.x = element_blank())
p1 | p2
We removed 692 (1.37%) outlier trials (150 ms < RT < 3000 ms).
df <- filter(df, Outlier == FALSE)
dfsub <- df |>
group_by(Participant) |>
select(Participant, Age, Sex, Education, Nationality, Ethnicity, Duration, Break_Duration, Screen_Resolution, Screen_Refresh, Device_OS) |>
slice(1) |>
ungroup()
The final sample included 46 participants (Mean age = 26.7, SD = 7.7, range: [19, 60]; Sex: 39.1% females, 56.5% males, 4.3% other; Education: Master, 17.39%; Bachelor, 41.30%; High School, 39.13%; Other, 2.17%).
plot_distribution <- function(dfsub, what = "Age", title = what, subtitle = "", fill = "orange") {
dfsub |>
ggplot(aes_string(x = what)) +
geom_density(fill = fill) +
geom_vline(xintercept = mean(dfsub[[what]]), color = "red", linetype = "dashed") +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
ggtitle(title, subtitle = subtitle) +
theme_modern() +
theme(
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(face = "italic", hjust = 0.5),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.y = element_blank()
)
}
plot_waffle <- function(dfsub, what = "Nationality") {
ggwaffle::waffle_iron(dfsub, what) |>
# mutate(label = emojifont::fontawesome('fa-twitter')) |>
ggplot(aes(x, y, fill = group)) +
ggwaffle::geom_waffle() +
# geom_point() +
# geom_text(aes(label=label), family='fontawesome-webfont', size=4) +
coord_equal() +
ggtitle(what) +
labs(fill = "") +
theme_void() +
theme(plot.title = element_text(face = "bold", hjust = 0.5))
}
p1 <- plot_distribution(dfsub, "Age", fill = "#FF9800")
p2 <- plot_distribution(dfsub, "Duration", title = "Total Duration", subtitle = "in minutes", fill = "#F44336")
p3 <- plot_distribution(dfsub, "Break_Duration", title = "Break Duration", subtitle = "in minutes", fill = "#3F51B5")
p4 <- plot_waffle(dfsub, "Sex") +
scale_fill_manual(values = c("Male" = "#2196F3", "Female" = "#E91E63", "Other" = "#FF9800"))
p5 <- plot_waffle(dfsub, "Education") +
scale_fill_viridis_d()
p6 <- plot_waffle(dfsub, "Nationality") +
scale_fill_metro_d()
p7 <- plot_waffle(dfsub, "Ethnicity") +
scale_fill_manual(values = c("Latino" = "#FF5722", "Asian" = "#FF9800", "Caucasian" = "#2196F3", "African" = "#4CAF50", "Jewish" = "#9C27B0"))
p8 <- plot_waffle(dfsub, "Screen_Resolution") +
scale_fill_pizza_d()
p9 <- plot_waffle(dfsub, "Device_OS") +
scale_fill_bluebrown_d()
# p10 <- plot_waffle(dfsub, "Screen_Refresh") +
# scale_fill_viridis_d()
(p1 / p2 / p3) | (p4 / p5 / p6) | (p7 / p8 / p9)
The analysis of study focused on the probability of errors as the main outcome variable. For each illusion, we started by visualizing the average effect of task difficulty and illusion strength to gain some insight on the underlying generative model. Next, we tested the performance of various models differing in their specifications, such as: with or without a log-transform on the task difficulty, with or without a 2nd order polynomial term of the illusion strength, and with or without the illusion side (up vs. down or left vs. right) as an additional predictor. We then fitted the best performing model under a Bayesian model, and compared its predicted links with that of a General Additive Model (GAM), which has an increased ability of mapping underlying potentially non-linear relationships at the expense of model simplicity. Finally, we used these models to estimate different values of task difficulty that would lead to similar error rates (2.5, 5 and 25%) when the illusion strength is null.
plot_descriptive <- function(data, side="leftright") {
if(side == "leftright") {
x <- data[data$Error == 0 & data$Illusion_Side == 1, ]$Answer[1]
x <- tools::toTitleCase(gsub("arrow", "", x))
if(x == "Left") {
data$Answer <- ifelse(data$Illusion_Side == 1, "Left", "Right")
} else if(x == "Right") {
data$Answer <- ifelse(data$Illusion_Side == 1, "Right", "Left")
}
} else {
x <- data[data$Error == 0 & data$Illusion_Side == 1, ]$Answer[1]
x <- tools::toTitleCase(gsub("arrow", "", x))
if(x == "Up") {
data$Answer <- ifelse(data$Illusion_Side == 1, "Up", "Down")
} else if(x == "Down") {
data$Answer <- ifelse(data$Illusion_Side == 1, "Down", "Up")
}
data$Answer <- fct_rev(data$Answer)
}
dodge1 <- 0.1 * diff(range(data$Illusion_Difference))
dodge2 <- -0.1 * diff(range(data$Illusion_Strength))
colors <- colorRampPalette(c("#4CAF50", "#009688", "#00BCD4", "#2196F3", "#3F51B5", "#673AB7", "#9C27B0"))(length(unique(data$Illusion_Strength)))
p1 <- data |>
group_by(Illusion_Difference, Illusion_Strength, Answer) |>
summarize(Error = mean(Error)) |>
mutate(Illusion_Strength = as.factor(round(Illusion_Strength, 2))) |>
ggplot(aes(x = Illusion_Difference, y = Error)) +
geom_bar(aes(fill=Illusion_Strength), position = position_dodge(width=dodge1), stat="identity") +
geom_line(aes(color = Illusion_Strength), position = position_dodge(width=dodge1)) +
scale_y_continuous(limits = c(0, 1), expand = c(0, 0), labels = scales::percent) +
scale_x_continuous(expand = c(0, 0)) +
scale_color_manual(values = colors) +
scale_fill_manual(values = colors) +
theme_modern() +
labs(
color = "Illusion Strength",
fill = "Illusion Strength",
y = "Probability of Error",
x = "Task Difficulty"
) +
theme(plot.title = element_text(face = "bold", hjust = 0.5))
colors <- colorRampPalette(c("#F44336", "#FFC107", "#4CAF50"))(length(unique(data$Illusion_Difference)))
p2 <- data |>
group_by(Illusion_Difference, Illusion_Strength, Answer) |>
summarize(Error = mean(Error)) |>
mutate(Illusion_Difference = as.factor(round(Illusion_Difference, 2))) |>
ggplot(aes(x = Illusion_Strength, y = Error)) +
geom_vline(xintercept=0, linetype="dotted", alpha=0.6) +
geom_bar(aes(fill=Illusion_Difference), position = position_dodge(width=dodge2), stat="identity") +
geom_line(aes(color = Illusion_Difference), position = position_dodge(width=dodge2)) +
scale_y_continuous(limits = c(0, 1), expand = c(0, 0), labels = scales::percent) +
scale_x_continuous(expand = c(0, 0)) +
scale_color_manual(values = colors) +
scale_fill_manual(values = colors) +
theme_modern() +
labs(
color = "Task Difficulty",
fill = "Task Difficulty",
y = "Probability of Error",
x = "Illusion Strength"
) +
theme(plot.title = element_text(face = "bold", hjust = 0.5))
if(side == "leftright") {
p <- ((p1 + facet_wrap(~Answer, ncol=2, labeller = "label_both")) /
(p2 + facet_wrap(~Answer, ncol=2, labeller = "label_both"))) +
plot_annotation(title = paste(data$Illusion_Type[1], "Illusion"),
theme = theme(plot.title = element_text(face = "bold", hjust = 0.5)))
} else {
p <- ((p1 + facet_wrap(~Answer, nrow=2, labeller = "label_both")) |
(p2 + facet_wrap(~Answer, nrow=2, labeller = "label_both"))) +
plot_annotation(title = paste(data$Illusion_Type[1], "Illusion"),
theme = theme(plot.title = element_text(face = "bold", hjust = 0.5)))
}
p
}
data <- filter(df, Illusion_Type == "Delboeuf")
plot_descriptive(data)
best_models <- function(data) {
models <- list()
for(i in 1:1) {
for(j in 1:1) {
for(k1 in c("", "_log", "_sqrt", "_cbrt")) {
for(k2 in c("")) {
for(side in c("", "-side")) {
name <- paste0("dif", k1, i, "-", "str", k2, j, side)
# print(name)
f <- paste0("poly(Illusion_Difference",
k1,
", ",
i,
") * poly(Illusion_Strength",
k2,
", ",
j,
") + (1|Participant)")
if(side == "-side") f <- paste0("Illusion_Side * ", f)
m <- glmmTMB::glmmTMB(as.formula(paste0("Error ~ ", f)),
data = data, family = "binomial")
if(performance::check_convergence(m)) {
models[[name]] <- m
}
}
}
}
}
}
to_keep <- compare_performance(models, metrics = c("BIC")) |>
arrange(BIC) |>
slice(1:10) |>
pull(Name)
test <- test_performance(models[to_keep], reference=1)
perf <- compare_performance(models[to_keep], metrics = c("BIC", "R2"))
merge(perf, test) |>
arrange(BIC) |>
select(Name, BIC, R2_marginal, BF) |>
mutate(BF = insight::format_bf(BF, name=""))
}
best_models(data)
## Name BIC R2_marginal BF
## 1 dif_cbrt1-str1 3305 0.392
## 2 dif_sqrt1-str1 3307 0.400 = 0.353
## 3 dif_log1-str1 3317 0.412 = 0.003
## 4 dif_cbrt1-str1-side 3330 0.394 < 0.001
## 5 dif_sqrt1-str1-side 3332 0.401 < 0.001
## 6 dif1-str1 3333 0.424 < 0.001
## 7 dif_log1-str1-side 3341 0.414 < 0.001
## 8 dif1-str1-side 3356 0.426 < 0.001
cbrt <- function(x) sign(x) * abs(x)**(1/3)
formula <- brms::bf(
Error ~ cbrt(Illusion_Difference) * Illusion_Strength +
(1 | Participant),
family = "bernoulli"
)
model <- brms::brm(formula,
data = data,
refresh = 0
)
## Finished in 0.6 seconds.
# parameters::parameters(model)
plot_model <- function(data, model) {
data <- mutate(data, .dots_side = ifelse(Error == 1, "bottom", "top"))
# Get variables
vars <- insight::find_predictors(model)$conditional
vardiff <- vars[1]
varstrength <- vars[2]
# Get predicted
pred <- estimate_relation(model,
at = vars,
length = c(NA, 25))
pred[[vardiff]] <- as.factor(pred[[vardiff]])
# Set colors for lines
colors <- colorRampPalette(c("#F44336", "#FFC107", "#4CAF50"))(length(unique(data[[vardiff]])))
diffvals <- as.numeric(as.character(unique(sort(pred[[vardiff]]))))
names(colors) <- diffvals
# Assign color from the same palette to every observation of data (for geom_dots)
closest <- diffvals[max.col(-abs(outer(data[[vars[1]]], diffvals, "-")))]
data$color <- colors[as.character(closest)]
data$color <- fct_reorder(data$color, closest)
# Manual jittering
xrange <- 0.05*diff(range(data[[varstrength]]))
data$x <- data[[varstrength]]
data$x[data$x > 0] <- data$x[data$x > 0] - runif(sum(data$x > 0), 0, xrange)
data$x[data$x < 0] <- data$x[data$x < 0] + runif(sum(data$x < 0), 0, xrange)
data$x[round(data$x, 2) == 0] <- data$x[round(data$x, 2) == 0] + runif(sum(round(data$x, 2) == 0), -xrange/2, xrange/2)
pred |>
ggplot(aes_string(x = varstrength, y = "Predicted")) +
geom_dots(
data = data,
aes(x=x, y = Error, group = Error, side = .dots_side, order=color),
fill = data$color,
color = NA,
alpha=0.5) +
geom_slab(data=data, aes(y = Error)) +
geom_ribbon(aes_string(ymin = "CI_low", ymax = "CI_high", fill = vardiff, group = vardiff), alpha = 0.2) +
geom_vline(xintercept = 0, linetype = "dashed") +
geom_hline(yintercept = c(0.05, 0.5, 0.95), linetype = "dotted", alpha = 0.5) +
geom_line(aes_string(color = vardiff, group = vardiff)) +
scale_y_continuous(limits = c(0, 1), expand = c(0, 0), labels = scales::percent) +
scale_x_continuous(expand = c(0, 0)) +
scale_color_manual(values = colors) +
scale_fill_manual(values = colors) +
coord_cartesian(xlim=c(min(data[[varstrength]]), max(data[[varstrength]]))) +
theme_modern() +
labs(
title = paste0(data$Illusion_Type[1], " Illusion"),
color = "Difficulty", fill = "Difficulty",
y = "Probability of Error",
x = "Illusion Strength"
) +
theme(plot.title = element_text(face = "bold", hjust = 0.5))
}
plot_model(data, model)
make_gam <- function(data) {
formula <- brms::bf(
Error ~ t2(Illusion_Difference, Illusion_Strength, bs = "cr", k=4) +
(1 | Participant),
family = "bernoulli"
)
model <- brms::brm(formula,
data = data,
refresh = 0
)
list(p = plot_model(data, model), model = model)
}
gam <- make_gam(data)
## Finished in 1.6 seconds.
gam$p
std_params <- function(model, min=0, max=2) {
estimate_relation(
model,
at = list(Illusion_Strength = c(0),
Illusion_Difference = seq(min, max, length.out=500)),
) |>
select(Illusion_Strength, Illusion_Difference, Error = Predicted) |>
slice(c(which.min(abs(Error - 0.005)),
which.min(abs(Error - 0.025)),
which.min(abs(Error - 0.25)))) |>
mutate(Error = insight::format_value(Error, as_percent=TRUE))
}
# range(data$Illusion_Difference)
# range(data$Illusion_Strength)
std_params(model, min=0, max=2)
## Illusion_Strength | Illusion_Difference | Error
## ------------------------------------------------
## 0.00 | 1.64 | 0.50%
## 0.00 | 0.77 | 2.51%
## 0.00 | 0.13 | 24.79%
std_params(gam$model, min=0, max=2)
## Illusion_Strength | Illusion_Difference | Error
## ------------------------------------------------
## 0.00 | 1.25 | 0.50%
## 0.00 | 0.68 | 2.50%
## 0.00 | 8.02e-03 | 25.01%
data <- filter(df, Illusion_Type == "Ebbinghaus")
plot_descriptive(data)
best_models(data)
## Name BIC R2_marginal BF
## 1 dif_cbrt1-str1 2212 0.615
## 2 dif_sqrt1-str1 2215 0.615 = 0.199
## 3 dif_cbrt1-str1-side 2227 0.623 < 0.001
## 4 dif_sqrt1-str1-side 2230 0.623 < 0.001
## 5 dif_log1-str1 2232 0.617 < 0.001
## 6 dif1-str1 2258 0.619 < 0.001
formula <- brms::bf(
Error ~ cbrt(Illusion_Difference) * Illusion_Strength +
(1 | Participant),
family = "bernoulli"
)
model <- brms::brm(formula,
data = data,
refresh = 0
)
## Finished in 0.7 seconds.
# parameters::parameters(model)
plot_model(data, model)
gam <- make_gam(data)
## Finished in 1.3 seconds.
gam$p
# range(data$Illusion_Difference)
# range(data$Illusion_Strength)
std_params(model, min=0, max=2)
## Illusion_Strength | Illusion_Difference | Error
## ------------------------------------------------
## 0.00 | 0.87 | 0.50%
## 0.00 | 0.44 | 2.48%
## 0.00 | 0.10 | 24.92%
std_params(gam$model, min=0, max=2)
## Illusion_Strength | Illusion_Difference | Error
## ------------------------------------------------
## 0.00 | 0.48 | 0.53%
## 0.00 | 0.44 | 2.85%
## 0.00 | 0.42 | 30.73%
data <- filter(df, Illusion_Type == "Rod-Frame")
plot_descriptive(data)
data <- filter(data, abs(Illusion_Strength) < 15)
plot_descriptive(data)
best_models(data)
## Name BIC R2_marginal BF
## 1 dif_sqrt1-str1 2070 0.464
## 2 dif_log1-str1 2072 0.458 = 0.371
## 3 dif_cbrt1-str1 2075 0.452 = 0.063
## 4 dif_sqrt1-str1-side 2095 0.478 < 0.001
## 5 dif1-str1 2095 0.505 < 0.001
## 6 dif_log1-str1-side 2097 0.471 < 0.001
## 7 dif_cbrt1-str1-side 2101 0.465 < 0.001
## 8 dif1-str1-side 2119 0.525 < 0.001
formula <- brms::bf(
Error ~ poly(Illusion_Difference, 2) * poly(Illusion_Strength, 2) +
(1 | Participant),
family = "bernoulli"
)
model <- brms::brm(formula,
data = data,
refresh = 0
)
## Finished in 0.7 seconds.
# parameters::parameters(model)
plot_model(data, model)
gam <- make_gam(data)
## Finished in 0.7 seconds.
gam$p
# range(data$Illusion_Difference)
# range(data$Illusion_Strength)
std_params(model, min=0.01, max=12)
## Illusion_Strength | Illusion_Difference | Error
## ------------------------------------------------
## 0.00 | 6.71 | 4.26%
## 0.00 | 6.71 | 4.26%
## 0.00 | 0.47 | 24.91%
std_params(gam$model, min=0.01, max=12)
## Illusion_Strength | Illusion_Difference | Error
## ------------------------------------------------
## 0.00 | 12.00 | 0.52%
## 0.00 | 4.53 | 2.51%
## 0.00 | 0.68 | 24.85%
data <- filter(df, Illusion_Type == "Vertical-Horizontal")
plot_descriptive(data)
data <- filter(data, abs(Illusion_Strength) < 90)
plot_descriptive(data)
best_models(data)
## Name BIC R2_marginal BF
## 1 dif_sqrt1-str1 2011 0.661
## 2 dif_cbrt1-str1 2014 0.657 = 0.234
## 3 dif_log1-str1 2017 0.674 = 0.050
## 4 dif1-str1 2022 0.679 = 0.004
## 5 dif_sqrt1-str1-side 2031 0.665 < 0.001
## 6 dif_cbrt1-str1-side 2034 0.662 < 0.001
## 7 dif_log1-str1-side 2038 0.677 < 0.001
## 8 dif1-str1-side 2043 0.681 < 0.001
formula <- brms::bf(
Error ~ sqrt(Illusion_Difference) * Illusion_Strength +
(1 | Participant),
family = "bernoulli"
)
model <- brms::brm(formula,
data = data,
refresh = 0
)
## Finished in 0.6 seconds.
# parameters::parameters(model)
plot_model(data, model)
gam <- make_gam(data)
## Finished in 0.6 seconds.
gam$p
# range(data$Illusion_Difference)
# range(data$Illusion_Strength)
std_params(model, min=0, max=0.40)
## Illusion_Strength | Illusion_Difference | Error
## ------------------------------------------------
## 0.00 | 0.32 | 0.50%
## 0.00 | 0.18 | 2.51%
## 0.00 | 0.04 | 24.93%
std_params(gam$model, min=0, max=0.40)
## Illusion_Strength | Illusion_Difference | Error
## ------------------------------------------------
## 0.00 | 0.30 | 0.50%
## 0.00 | 0.20 | 2.50%
## 0.00 | 0.04 | 24.89%
data <- filter(df, Illusion_Type == "Zöllner")
plot_descriptive(data)
data <- filter(data, abs(Illusion_Strength) < 45)
plot_descriptive(data)
best_models(data)
## Name BIC R2_marginal BF
## 1 dif_cbrt1-str1 1041 0.405
## 2 dif_log1-str1 1043 0.414 = 0.464
## 3 dif_sqrt1-str1 1044 0.423 = 0.260
## 4 dif1-str1 1065 0.491 < 0.001
## 5 dif_cbrt1-str1-side 1066 0.409 < 0.001
## 6 dif_log1-str1-side 1068 0.418 < 0.001
## 7 dif_sqrt1-str1-side 1069 0.427 < 0.001
## 8 dif1-str1-side 1090 0.494 < 0.001
formula <- brms::bf(
Error ~ poly(Illusion_Difference, 2) * Illusion_Strength +
(1 | Participant),
family = "bernoulli"
)
model <- brms::brm(formula,
data = data,
refresh = 0
)
## Finished in 0.4 seconds.
# parameters::parameters(model)
plot_model(data, model)
# range(data$Illusion_Difference)
# range(data$Illusion_Strength)
std_params(model, min=0, max=5)
## Illusion_Strength | Illusion_Difference | Error
## ------------------------------------------------
## 0.00 | 3.55 | 3.04%
## 0.00 | 3.55 | 3.04%
## 0.00 | 0.22 | 25.07%
data <- filter(df, Illusion_Type == "White")
plot_descriptive(data)
best_models(data)
## Name BIC R2_marginal BF
## 1 dif_log1-str1 2176 0.783
## 2 dif_cbrt1-str1 2178 0.784 = 0.404
## 3 dif_sqrt1-str1 2181 0.784 = 0.101
## 4 dif_log1-str1-side 2188 0.788 = 0.003
## 5 dif_cbrt1-str1-side 2190 0.789 = 0.001
## 6 dif_sqrt1-str1-side 2192 0.790 < 0.001
## 7 dif1-str1 2195 0.787 < 0.001
## 8 dif1-str1-side 2206 0.793 < 0.001
formula <- brms::bf(
Error ~ log(Illusion_Difference) * poly(Illusion_Strength, 2) +
(1 | Participant),
family = "bernoulli"
)
model <- brms::brm(formula,
data = data,
refresh = 0
)
## Finished in 0.8 seconds.
# parameters::parameters(model)
plot_model(data, model)
gam <- make_gam(data)
## Finished in 0.8 seconds.
gam$p
# range(data$Illusion_Difference)
# range(data$Illusion_Strength)
# unique(data$Illusion_Strength)
std_params(model, min=0, max=20)
## Illusion_Strength | Illusion_Difference | Error
## ------------------------------------------------
## 0.00 | 20.00 | 0.71%
## 0.00 | 12.75 | 2.50%
## 0.00 | 5.13 | 25.10%
std_params(gam$model, min=0, max=20)
## Illusion_Strength | Illusion_Difference | Error
## ------------------------------------------------
## 0.00 | 9.06 | 0.50%
## 0.00 | 5.45 | 2.50%
## 0.00 | 0.84 | 25.06%
data <- filter(df, Illusion_Type == "Müller-Lyer")
plot_descriptive(data, side = "updown")
best_models(data)
## Random effect variances not available. Returned R2 does not account for random effects.
## Random effect variances not available. Returned R2 does not account for random effects.
## Name BIC R2_marginal BF
## 1 dif_cbrt1-str1 2449 0.749
## 2 dif_log1-str1-side 2534 0.748 < 0.001
## 3 dif_cbrt1-str1-side 2535 0.757 < 0.001
formula <- brms::bf(
Error ~ log(Illusion_Difference) * poly(Illusion_Strength, 2) +
(1 | Participant),
family = "bernoulli"
)
model <- brms::brm(formula,
data = data,
refresh = 0
)
## Finished in 0.8 seconds.
# parameters::parameters(model)
plot_model(data, model)
gam <- make_gam(data)
## Finished in 1.0 seconds.
gam$p
# range(data$Illusion_Difference)
# range(data$Illusion_Strength)
std_params(model, min=0, max=0.6)
## Illusion_Strength | Illusion_Difference | Error
## ------------------------------------------------
## 0.00 | 0.60 | 2.66%
## 0.00 | 0.60 | 2.66%
## 0.00 | 0.07 | 24.86%
std_params(gam$model, min=0, max=0.6)
## Illusion_Strength | Illusion_Difference | Error
## ------------------------------------------------
## 0.00 | 0.45 | 0.50%
## 0.00 | 0.27 | 2.51%
## 0.00 | 0.02 | 24.89%
data <- filter(df, Illusion_Type == "Ponzo")
plot_descriptive(data, side = "updown")
best_models(data)
## Name BIC R2_marginal BF
## 1 dif_cbrt1-str1 2551 0.661
## 2 dif_sqrt1-str1 2555 0.665 = 0.187
## 3 dif_log1-str1 2576 0.678 < 0.001
## 4 dif1-str1 2590 0.684 < 0.001
## 5 dif_log1-str1-side 2592 0.687 < 0.001
## 6 dif1-str1-side 2606 0.694 < 0.001
formula <- brms::bf(
Error ~ poly(Illusion_Difference, 2) * poly(Illusion_Strength, 2) +
(1 | Participant),
family = "bernoulli"
)
model <- brms::brm(formula,
data = data,
refresh = 0
)
## Finished in 0.9 seconds.
# parameters::parameters(model)
plot_model(data, model)
gam <- make_gam(data)
## Finished in 0.8 seconds.
gam$p
# range(data$Illusion_Difference)
# range(data$Illusion_Strength)
std_params(model, min=0, max=20)
## Illusion_Strength | Illusion_Difference | Error
## ------------------------------------------------
## 0.00 | 0.36 | 3.17%
## 0.00 | 0.36 | 3.17%
## 0.00 | 0.04 | 22.41%
std_params(gam$model, min=0, max=20)
## Illusion_Strength | Illusion_Difference | Error
## ------------------------------------------------
## 0.00 | 9.14 | 0.50%
## 0.00 | 0.24 | 2.15%
## 0.00 | 0.04 | 30.50%
data <- filter(df, Illusion_Type == "Poggendorff")
plot_descriptive(data, side = "updown")
data <- filter(data, abs(Illusion_Strength) < 45)
plot_descriptive(data, side = "updown")
best_models(data)
## Name BIC R2_marginal BF
## 1 dif_cbrt1-str1 1900 0.527
## 2 dif_sqrt1-str1 1905 0.526 = 0.086
## 3 dif_log1-str1 1923 0.522 < 0.001
## 4 dif_cbrt1-str1-side 1929 0.529 < 0.001
## 5 dif1-str1 1930 0.521 < 0.001
## 6 dif_sqrt1-str1-side 1934 0.529 < 0.001
## 7 dif_log1-str1-side 1952 0.527 < 0.001
## 8 dif1-str1-side 1958 0.526 < 0.001
formula <- brms::bf(
Error ~ log(Illusion_Difference) * poly(Illusion_Strength, 2) +
(1 | Participant),
family = "bernoulli"
)
model <- brms::brm(formula,
data = data,
refresh = 0
)
## Finished in 0.5 seconds.
# parameters::parameters(model)
plot_model(data, model)
gam <- make_gam(data)
## Finished in 0.6 seconds.
gam$p
# range(data$Illusion_Difference)
# range(data$Illusion_Strength)
std_params(model, min=0, max=0.5)
## Illusion_Strength | Illusion_Difference | Error
## ------------------------------------------------
## 0.00 | 0.50 | 1.20%
## 0.00 | 0.21 | 2.50%
## 0.00 | 0.01 | 25.40%
std_params(gam$model, min=0, max=0.5)
## Illusion_Strength | Illusion_Difference | Error
## ------------------------------------------------
## 0.00 | 0.28 | 1.65%
## 0.00 | 0.16 | 2.50%
## 0.00 | 0.00 | 12.19%
data <- filter(df, Illusion_Type == "Contrast")
plot_descriptive(data, side = "updown")
best_models(data)
## Random effect variances not available. Returned R2 does not account for random effects.
## Random effect variances not available. Returned R2 does not account for random effects.
## Name BIC R2_marginal BF
## 1 dif_sqrt1-str1 2353 0.673
## 2 dif_cbrt1-str1 2354 0.673 = 0.713
## 3 dif_log1-str1 2357 0.672 = 0.128
## 4 dif1-str1 2359 0.675 = 0.047
## 5 dif_cbrt1-str1-side 2580 0.685 < 0.001
## 6 dif1-str1-side 2585 0.685 < 0.001
formula <- brms::bf(
Error ~ Illusion_Difference * poly(Illusion_Strength, 2) +
(1 | Participant),
family = "bernoulli"
)
model <- brms::brm(formula,
data = data,
refresh = 0
)
## Finished in 0.6 seconds.
plot_model(data, model)
gam <- make_gam(data)
## Finished in 0.8 seconds.
gam$p
# range(data$Illusion_Difference)
# range(data$Illusion_Strength)
std_params(model, min=0, max=25)
## Illusion_Strength | Illusion_Difference | Error
## ------------------------------------------------
## 0.00 | 15.48 | 0.50%
## 0.00 | 11.72 | 2.50%
## 0.00 | 5.81 | 24.87%
std_params(gam$model, min=0, max=25)
## Illusion_Strength | Illusion_Difference | Error
## -----------------------------------------------
## 0.00 | 18.79 | 0.50%
## 0.00 | 9.32 | 2.51%
## 0.00 | 0.00 | 9.54%
This study allowed to elaborate a more precise prior idea regarding the magnitude of the effects at stake and the type of interaction between them. Crucially, it allowed us to refine the range of task difficulty and illusion strength for the stimuli of the next study to maximize the information gain.